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ABSTRACT 

In  this  paper  we  will  develop  a  method  to  determine 
cross  sections  of  arbitrary  two-dimensional  tubular  struc¬ 
tures,  which  are  allowed  to  branch,  by  means  of  a  Stokes 
flow  based  boundary  integral  formulation.  The  measure  for 
the  cross  sections  for  a  point  on  the  boundary  of  a  given 
structure  will  be  the  path  obtained  by  integrating  perpen¬ 
dicularly  to  the  flow  lines  from  one  side  of  the  boundary  to 
the  other.  Special  emphasis  will  be  put  on  the  behavior  at 
branching  points,  the  behavior  at  vortices,  and  the  necessary 
boundary  conditions.  The  method  can  be  extended  to  three 
dimensional  problems. 

1.  INTRODUCTION 

Measuring  cross  sections  of  structures  in  a  consistent, 
meaningful  way  is  important  for  many  applications.  In 
medical  imaging,  anatomical  structures  often  times  exhibit 
complicated,  convoluted  shapes.  Manual  thickness  mea¬ 
surements  based  on  image  data  (especially  in  three  dimen¬ 
sions)  then  easily  become  error-prone.  Jones  et  al.  [1]  use 
Laplace’s  equation  to  measure  cortical  thickness  in  three- 
dimensional  images  whose  variations  can  be  associated 
with  many  pathologies:  e.g.  Alzheimer’s  disease.  Yezzi 
and  Prince  [2]  improve  the  speed  of  Jones’  computational 
method  and  eliminate  the  need  for  the  explicit  computation 
of  trajectories.  However,  as  in  [1],  the  approach  is  confined 
to  structures  that  can  be  described  by  two  simply  connected 
boundaries.  Budding  structures  (see  Section  5)  pose  se¬ 
rious  problems  for  the  aforementioned  approaches.  Note 
that,  since  the  temperature  at  the  two  boundaries  is  fixed, 
no  sensible  thickness  would  be  assigned  to  the  cavity  of  the 
budding  structure  by  a  Laplace  equation  based  method:  a 
solution  of  Laplace’s  equation  will  not  exhibit  vortices. 

This  paper  will  deal  with  structures  that  cannot  be  de¬ 
scribed  by  two  simply  connected  boundaries,  e.g.  struc¬ 
tures  that  branch  (e.g.  blood  vessels),  and/or  exhibit  bud¬ 
ding/constriction  of  boundary  part(s).  Instead  of  defining 
thickness  by  means  of  the  solution  of  Laplace’s  equation 
we  use  the  Stokes  equation  (describing  fluid  flow  at  low 
Reynolds  numbers),  with  free  slip  boundary  conditions.  We 
assume  that  the  object  to  be  measured  is  given  as  a  digital 
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image.  To  solve  the  Stokes  equation  we  use  a  boundary  in¬ 
tegral  formulation,  thus  reducing  the  dimensionality  of  the 
problem  by  one  dimension.  There  is  no  need  to  discretize 
the  interior  of  the  object.  We  restrict  ourselves  to  the  two  di¬ 
mensional  problem  (flow  in  a  plane).  However,  the  method¬ 
ology  developed  can  be  extended  to  three  dimensions  (the 
boundary  integral  formulation  for  the  Stokes  equation  is 
standard  in  three  dimensions).  Other  fluid  flow  equations 
have  recently  been  used  in  computer  graphics  and  image 
processing,  notably  for  image  inpainting/reconstruction  [3], 
Examples  for  using  the  boundary  element  approach  in  com¬ 
puter  vision  can  be  found  in  [4], 

Section  2  introduces  the  Stokes  equations  and 
presents  the  boundary  integral  formulation  as  given 
by  Pozrikidis  [5],  Section  3  deals  with  the  discretization 
of  the  boundary  integral  equations,  and  introduces  the 
formalism  for  solving  the  discretized  equations.  Section  4 
describes  the  methodology  for  measuring  cross  sections 
once  a  flow  field  has  been  obtained.  Examples  are  given 
in  Section  5.  The  paper  ends  with  a  conclusion  and 
suggestions  for  future  work  in  Section  6. 

2.  BOUNDARY  INTEGRAL  FORMULATION  OF 
THE  STOKES  EQUATIONS 

Pozrikidis  [5]  gives  a  good  introduction  to  the  boundary 
integral  method  for  linearized  viscous  flow.  Higdon  [6] 
treats  two-dimensional  Stokes  flow  problems  for  arbitrarily 
shaped  domains.  Zeb  et  al.  [7]  present  an  implementation  of 
Higdon’s  approach  and  extend  it  to  handle  pressure  bound¬ 
ary  conditions.  This  section  aims  at  reviewing  the  basics 
for  the  formulation  of  the  Stokes  flow  problem  in  terms  of 
boundary  integrals,  following  [5,  7]. 

2.1.  The  Stokes  flow  equations 

Incompressible  fluid  flow  problems  where  viscous  forces 
dominate  can  be  described  by  the  divergence-free  condition 
and  the  Stokes  equation 

V  •  u  =  0,  —  VP  +  [iX72u  +  pb  =  0.  (1) 

Here  u  is  the  fluid  velocity,  //  the  viscosity  of  the  fluid,  P  is 
the  pressure  and  b  is  the  body  force.  V  indicates  the  gra¬ 
dient,  V-  represents  the  divergence  operator,  and  V2  the 
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Laplacian.  In  the  sequel  we  incorporate  the  body  force  in 
the  modified  pressure  pMOD  p  pb  ■  x.  To  avoid  un¬ 
necessary  heavy  notation  P  is  understood  to  mean  PMOD 
in  what  follows. 


2.2.  Boundary  integral  formulation 


The  incompressible  solution  of  the  Stokes  equation 
—VP  +  gS/2u  +  gS  (x  —  xq)  =  0, 


where  5(x  —  Xq)  is  a  two-dimensional  Dirac-delta-function 
centered  at  Xo  and  g  is  an  arbitrary  vector-valued  constant, 
is  given  by1 


Ui(x) 


47 TJU 


Gij  (a?,  Xq 'jgj . 


(2) 


(2)  is  the  solution  for  the  Stokes  flow  caused  by  a  two- 
dimensional  point  force  at  Xq.  Gij  are  the  Green’s  func¬ 
tions,  which  describe  the  influence  of  a  force  at  position  Xq 
on  the  solution  at  position  x.  Since  the  Stokes  flow  prob¬ 
lem  is  linear,  one  can  express  the  solution  to  an  arbitrary 
force  distribution  by  superposition  of  the  respective  Green’s 
function  solutions  for  each  point  force. 

The  stress  associated  with  this  flow  may  be  written 


&ij  (xj  I'ij  k(x .  XQ^(jk. 
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where  the  stress  tensor  Tijk  is 


Tijki^X^X  o)  — 

£  (  A  dGij  (  ^  dGkj  (  A 

-  6ikPj{X,  Xq)  +  -Q^p-yx,  Xq)  +  -£-±(x,X0). 

The  Green’s  functions  Gij  (also  called  the  two-dimensional 
Stokeslets)  and  the  stress  tensor  Tijk  are  given  by  [5] 


G^  =  -dij  in  r 
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3.  BOUNDARY  DISCRETIZATION 

Following  [7]  we  assume  a  piecewise  linear  parameteriza¬ 
tion  of  the  boundary.  The  boundary  dD  is  divided  into  n 
straight  lines  dDk ,  k  =  1,2, . .  .  n,  where  fi(x)  and  Ui(x) 
are  assumed  constant  along  these  elements.  With  these  sim¬ 
plifications,  the  governing  boundary  integral  equations  (3) 
become 
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where  the  two  dimensional  Stokeslet  integrals  are  given  by 
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/  —  in  r  + 
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and  the  stress  tensor  integrals  by 


(5) 


Axjn , 

f 

X\X\ 

di 

ldDk 

r * 

Aiin , 

i(Xk) 

f 

X\X2 

d  l  = 

=  T21k 

ldDk 

r * 

Axin, 

i(Xk) 

f 

X2X2 

d  l. 

(6) 

ldDk 

r * 

Here  6lk  is  the  Kronecker-delta-function,  r  =  \\x  —  atoll 2, 
and  x  =  x  —  Xq. 

For  a  point  xq  inside  or  on  the  boundary  the  solution  to 
the  Stokes  equation  (1)  in  terms  of  boundary  integral  equa¬ 
tions  is 

r](xo)uj(xo)  =  f  fi(x)Gij(x,x 0)  d l(x)  + 

dTTft  J  dD 

/  ,Uji(x')Tijk(x ,  d /(at),  (3) 
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for  the  velocities,  where  fi(x)  =  aij(x)n,j(x)  is  the  stress 
force,  and  n(x)  is  the  normal  vector  to  the  boundary  OD 
at  position  x,  pointing  towards  the  interior  of  the  solution 
domain  D.  tj(xo)  accounts  for  the  dependency  of  the  ve¬ 
locity  integral  on  the  location  of  Xq:  ij(xq)  =  1  if  Xq  £  D, 
g(x0)  =  \  if  x0  g  dD. 

1  Note  that  we  are  looking  at  the  two-dimensional  Stokes  flow  problem. 
This  is  the  solution  of  the  two-dimensional  equation,  written  with  Ein¬ 
stein’s  summation  convention  (summing  over  multiply  occurring  indices). 
Solutions  in  three  dimensions  also  exist,  as  treated  in  [5].  We  denote  vector 
components  by  subscripts,  i.e.  u\  and  U2  for  the  components  of  u. 


(\I  is  the  differential  along  the  boundary  element.  Equa¬ 
tions  (5,  6)  can  be  evaluated  analytically;  they  describe  the 
influence  of  a  whole  boundary  segment  on  a  point  ato  lo¬ 
cated  in  or  on  the  boundary  of  the  domain  D.  The  inte¬ 
grands  of  Equations  (5,  6)  exhibit  singularities.  These  are 
removable  (except  for  the  case  of  points  which  lie  exactly 
on  the  end  point  of  a  boundary  element)2. 

To  evaluate  the  discretized  boundary  integrals  (4)  we  use 
the  following  two  step  approach  (see  [7]): 

a)  Evaluate  the  velocity  boundary  integrals  for  the  n 
center  points  of  the  boundary  elements,  by  specify¬ 
ing  2 n  boundary  conditions.  This  is  a  linear  equa¬ 
tion  system  with  An  unknowns  (/i,  /2,  Ui,  w2  for  ev¬ 
ery  boundary  element)  and  2 n  equations.  Fast  algo¬ 
rithms  for  solving  such  boundary  integral  equations 
exist,  e.g.  [8], 

b)  Solve  for  the  velocity  at  arbitrary  positions  inside  the 
domain  D  by  using  Equation  (4). 

2We  omit  the  explicit  analytical  expressions  for  the  integrals  due  to 
space  constraints. 


The  boundary  decomposes  into  the  inlet,  outlet  and  wall 
boundary  parts:  dD  =  dI),\JdI)0lJ()Dn!.  We  prescribe  the 
velocity  profile  ui,  U2 3  at  the  inlet  ( dDi )  and  assume  stress 
free  boundary  conditions  for  the  outlet(s)  (i.e.  /i  =  J  2  =  0 
for  x  £  dD,,).  dDw  exhibits  perfect  slip  without  perme¬ 
ation:  U2  =  0,  /  j  =  0.  Note,  that  these  boundary  condi¬ 
tions  directly  translate  to  the  three-dimensional  case.  While 
the  Stokes  equation  itself  is  linear  and  reversible,  we  note 
that  the  generated  flow  fields  will  vary  slightly  under  inter¬ 
changing  the  identification  of  the  inlet  and  outlet,  because 
of  the  prescribed  boundary  conditions.  These  differences, 
however,  will  largely  be  restricted  to  the  immediate  vicini¬ 
ties  of  the  inlet  and  outlet  for  reasonable  prescribed  inlet 
flows. 

4.  MEASUREMENT  OF  CROSS  SECTIONS 

To  measure  cross  sections  of  a  structure  we  propose  to  inte¬ 
grate  perpendicular  to  the  computed  flow  field.  We  compute 

x(l)  =  f  v  dl,  (7) 

Jo 

where  v  £  K2  is  a  unit  vector  normal  to  the  fluid  speed  u  at 
all  points,  x(L)  £  dD  with  I/O,  and  a;(0)  is  the  bound¬ 
ary  point  for  which  the  cross  section  is  to  be  measured.  L 
is  then  the  path  length  from  the  initial  to  the  final  boundary 
point.  We  set  the  initial  condition  of  the  vector  to  the  nor¬ 
mal  vector  for  the  boundary  element  of  the  starting  point, 
and  disallow  directional  flips  throughout  the  integration. 

For  the  evaluation  of  (7)  we  use  a  modified  version 
of  the  variable  order  Runge-Kutta  method  by  Cash  and 
Karp  [9].  It  is  modified  in  the  sense  that  a  minimal  step  size 
hmin  is  introduced.  Once  the  integration  algorithm  tries  to 
use  a  step  size  h  <  hmin  we  stop  the  integration  and  de¬ 
clare  the  current  integrator  state  as  the  end  point.  In  this 
way  we  stop  the  integration  when  either  a  boundary  point 
or  a  stagnant  point  (where  the  step  size  will  tend  to  zero) 
is  reached.  The  combination  of  an  integrator  with  adap¬ 
tive  step  size  and  the  boundary  integral  method  thus  yields 
a  method  that  will  naturally  adapt  to  the  problem  at  hand. 
If  a  lot  is  known  about  the  shape  of  the  boundary  (i.e.  we 
have  a  high  resolution  image)  this  information  will  be  used, 
without  requiring  too  many  steps  to  integrate  from  an  ini¬ 
tial  point  on  the  boundary  to  the  target  point.  For  very  thin 
structures  (thickness  in  the  order  of  a  few  pixels)  the  method 
will  automatically  result  in  subpixel  accuracy. 

5.  EXAMPLES 

This  section  presents  examples  for  cross  section  measure¬ 
ments:  a  budding  domain,  and  a  spinous  process  of  a  ver¬ 
tebra.  We  set  fi  =  1  and  hmin  =  0.005,  and  the  relative 
error  tolerance  of  the  integrator  to  e  =  10-4.  Figures  1(a) 
shows  the  discretized  domain  of  the  budding  example,  with 


(a)  Discretization,  with  boundary  elements 
and  associated  normal  vectors. 


(b)  Normalized  velocity  field.  Light  gray  ar¬ 
rows  indicate  the  prescribed  inflow  veloci¬ 
ties. 


(c)  Cross  sections. 

Fig.  1.  The  budding  domain. 


its  boundary  elements  and  the  associated  normal  vectors. 
Figures  1(b)  and  2(a)  show  the  normalized  flow  fields  of 
the  examples,  and  the  computer  tomography  (CT)  image  of 
the  spinous  process4.  Figures  1(c)  and  2(b)  show  their  re¬ 
spective  measures  for  the  cross  sections.  Of  special  interest 
is  the  behavior  close  to  the  branching  point  for  the  spinous 
process.  While  the  Stokes  flow  based  approach  produced 
sensible  results,  a  method  based  for  example  on  finding  the 
nearest  point  on  the  opposite  boundary  would  lead  to  unrea¬ 
sonable  results.  The  budding  domain  exhibits  a  vortex  in  its 
flow  field  (see  Figure  1(b)).  The  path  lines  perpendicular  to 
the  flow  field  of  the  budding  domain  capture  the  rectangular 
part  of  the  domain  well.  The  cavity  shows  the  end  points 
of  several  path  lines  at  the  stagnant  center  of  the  vortex. 


3  Overlined  variables  denote  variables  given  in  a  locally  defined  coordi¬ 
nate  system  of  a  boundary  element,  where  the  X2  direction  coincides  with 
the  normal  vector. 


4The  authors  would  like  to  thank  Prof.  Yamamoto  for  providing  the 
image.  The  boundaries  were  extracted  by  means  of  an  active  contours 
based  approach. 


6.  CONCLUSION  AND  FUTURE  WORK 


(a)  Normalized  velocity  field  and  CT  image.  Light 
gray  arrows  indicate  the  prescribed  inflow  veloci¬ 
ties. 


Fig.  2.  The  spinous  process. 


One  could  for  example  connect  these  path  lines  by  requir¬ 
ing  minimal  change  of  direction  at  the  joining  point.  We 
see  that  through  the  presence  of  the  vortex,  a  measure  for 
the  cross  section  of  the  budding  part  of  the  structure  as  well 
as  of  the  rectangular  part  can  be  given.  This  would  not  have 
been  possible  with  an  approach  based  on  Laplace’s  equation 
(heat  flow).  Note,  that  the  observations  made  in  this  section 
will  have  correspondences  in  the  three-dimensional  case. 


In  this  paper  a  Stokes  flow  based  method  for  measuring 
cross  sections  of  two-dimensional  structures  was  presented. 
This  method  is  more  versatile  than  the  Laplace  equation 
based  approach.  A  greater  variety  of  domains  can  be  han¬ 
dled;  e.g.  branching  and  budding  domains.  There  is  no  re¬ 
striction  to  domains  that  are  enclosed  by  two  simply  con¬ 
nected  surfaces.  A  boundary  integral  approach  was  used 
to  deal  easily  with  arbitrary  shaped  boundaries,  yielding  a 
clean  method,  without  the  need  for  discretization  of  the  in¬ 
terior  of  the  domain.  The  method  is  extendable  to  three 
dimensions,  where  one  of  the  natural  geometric  measure¬ 
ments  would  be  cross-sectional  surface  area,  and  we  will 
expect  to  find  lines  of  stagnant  points  (instead  of  isolated 
stagnant  points  at  the  center  of  vortices).  Future  work  could 
include  the  extension  to  three  dimensional  problems,  inves¬ 
tigations  on  the  handling  of  stagnant  points/lines  (the  vor¬ 
tices),  and  center  line  construction. 
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